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Abstract. First we present the general equation form of a thermal explo¬ 
sion in a vessel with boundary values, later use central difference method and 
Newton iteration method to solve the relevant partial differential equations in 
one-dimensional and two-dimensional forms, finally we use numerical experi¬ 
ments to verify the order of convergence of the central difference scheme and 
provide the experiment results. 


Introduction 


N.N.Semenov [1], Ya.B.Zeldovich [2], D.A.Frank-Kameneetiskiaei [3], O.M.Todes 
and P.V.Melent’ev [7 , A.G.Merzhanov and F.I.Dubovitsky |5j, B.Gray [B] had done 
a lot of research work in the field of thermal explosion theory. 

An example of a boundary value problem(BVP) describing the onset of a thermal 
explosion in a vessel is the following [2]: 


f aV 2 T + Qk(T) = 0,xefi; 
{ T(x) =T 0 ,x€ dfl 


where T(x) is the temperature in the vessel([T] = K). Equation ( |0.0.1 ) is a 
reaction-diffusion equation and determines the balance between conduction(thermal 
diffusion) and heat production by chemical reactions. Other variables/parameters 
in the BVP ( 0.0. 1| ) are the ambient temperature T 0 ([T 0 ] = K), the thermal diffu- 
sivity a ([a] = m 2 /s), the heat release parameter Q([Q] = K) and the reaction rate 
k = k(T)([k\ = s- 1 ), given by the expression 


k(T) = Ae~ TJT , 

where A is the pre-exponential factor ([A] = s _1 ) and T a the activation temperature 
([T a ] = K). It is customary to rewrite equation |0. 0.1 in terms of the dimensionless 
variables u and x*, defined by 

f T a T(x) - To * X 

u{x ) := —-—- ,x := —, 


Tn 


T„ 
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where £ is a characteristic dimension of the vessel. To first order approximation, 
We can easily obtain 


( 0 . 0 . 2 ) 


V 2 u - 


qe 


= 0 , 


9 = 


QA£ 2 T a e 


-e 


aTn 


T 

9 = —. 
T 0 


A solution of equation (0.0.2) exists if q is small enough, i.e.,g < q* for some 
critical value q* , meaning that chemical heat production can be entirely balanced 


by conduction. However, for increasing q the solution of equation (0.0.2) suddenly 


does not exist anymore. This phenomenon is referred to as a thermal explosion. 


Since (0.0.2) is a nonliear problem, [7] suggests we should use one of the following 
approaches: Newton iteration, Gauss-Jacobi iteration and transient methods. In 
this paper, we would first use central difference scheme to discretize the nonlinear 


equation (0.0.2), later use Newton iteration method 0191 HQj to solve the nonlinear 
system of equations. 

The paper is organized as the following: in section |TJ the one-dimensional case 
of equation (0.0.2) is proposed and the solution is given analytically, later we use 


central difference method and Newton iteration method as the numerical algorithm 
to solve this equation, specially we use numerical experiments to obtain the thresh¬ 
old value of parameter q and the convergence order of the algorithm by Richardson 
extrapolation method; in section [2j the two-dimensional case of equation (0.0.2) is 
proposed and again we use the same algorithm as in section [l] to solve the equation 
numerically with different values of parameter g, the numerical experiment results 
are listed as a table and several figures; finally the conclusion is given in section [3] 

1. Solve One-dimensional thermal Explosion Model 


Consider the following one-dimensional BVP of (0.0.2) 

n , / u" + qe u = 0,x G (0,1), 

1 1 u'(0) = 0,u(l)=0. 


In this section, we will give the analytic of equation (1.0.3), later use central 


difference method and Newton iteration method to compute the numerical solution 
of this equation, finally use Richardson extrapolation to compute the convergence 
order of the numerical scheme. 


1.1. Analytic solution. The solution of equation (1.0.3) can be found by subse¬ 


quently multiplying equation with u', integrate, isolate u' and integrate again. This 
solution is given by 


( 1 . 1 . 1 ) 


2 m 2 

u(x ) = ln{ -) — 2 ln(cosh(p,x)), 

9 


where the parameter p, satisfies the relation 

cosh/j, = \J 2 fq 11 . 

1.2. Central difference and Newton iteration Scheme. In order to compute 
the numerical solution of equation (|1.0.3 ), we introduce the grid 


Xj = (j - l)Aa;, j = 1,2, • • • , M, 

where Ax = 1/{M — 1) is the grid size. We first use standard central differences 
for all derivatives in equation (1.0.3). Let Uj denote the numerical approximation 
of u(xj) and let u := ('u 1 ,u 2 ,--- ,um- 1 ) be the vector of unknowns. Specially, 
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we introduce the virtual point xq and apply the central difference to discretise the 
Robin boundary condition at the left boundary as 


Thus we have 


u 2 - Up 

2Ax 


u'(0) = 0 


MO = U 2 - 

The resulting nonlinear system of algebraic equations can be written in the form 

(1.2.1) N(u) := Au + f(u) = 0, 

where 

(-2 2 0 • • • 0 \ 

1 -2 1 : 


f(u) = q ■ diag(e Ul , e“ 2 ,• • • ,e UM ~ 1 ),A = 


1 


Ax 2 


0 


\0 


0 


0 1 - 2 / 


Now we would use Newton iteration scheme to solve equation (1.2.1). First give a 
suitable initial guess. Solve 

i" + q = 0,2 <E (0,1), 


we obtain 


T(0) = 0, m(1) = 0, 




The Newton iteration scheme to solve equation (1.2.1) is the following: 


with 


given, u 0 , 

solve J(u l )s l = —N(u l ), 
u i+i = u i + s i ; 


u° = |(1 -x 2 ),J(u) 


dN(u) 

du 


, dN z (u) 

duj 


),s l 


N(u l ) 
J{u l )' 


1.3. Numerical scheme to compute threshold value of Parameter q. The 

numerical experiment result shows that the threshold value of parameter q is 0.878. 
The numerical scheme to compute threshold value of Parameter q is the following: 
Input:grid size Aa;,maxit(maximum iteration),tol(tolerance). 

Parameters:g = 0.87 : 0.001 : 0.88. 

Output:threshold value of parameter q*. 

M°:=f(l-z 2 ); 

for / =: 1 to length(g) do 

solve J(u l )s l = 
u i+i _ u i _|_ s i. 

F = A*u l+1 + f(u l+1 ); 
nF = norm{F, 2); 
con v=(nF < tol\(l == 
end 

if (nF > tol) 

q * = q(i - i); 


maxit))\ 
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break 

else 

continue 

end 

* 

Q 


1.4. Convergence order. Due to the nonlinearity of the problem, convergence 
is difficult to prove. Instead, we verify the order of convergence by numerical 
experiments. Verify the order of convergence of the central difference scheme by 


computing numerical solutions of equation (1.2.11 for several values of the grid size 
Ax. Using Richardson Extrapolation method, the order of convergence can be 
computed by 


ln{ M2 ~ Me °° ) 
' Wuo— ' 


P = 


||U2-U e 

ii «2- 

ln2 


where u e represents the analytic solution from equation (1.1.1), zti,u 2 represents 


numerical solutions from different grid sizes. Numerical experiments show that the 
order of convergence is 2. 

2. Solve One-dimensional thermal Explosion Model 

We consider the following two-dimensional BVP: 

fyr + + qe u = 0, x £ (0, £), y G (0,1), 

u{0,y) = 0 ,u(£,y) = g{y),y G (0,1), 
x G (0,£), 


( 2 . 0 . 1 ) 


^(x,0) = f(x,l) = 0, 


representing a rectangular vessel of aspect ratio i with the horizontal walls isolated. 
For the boundary function g(y) we take 

0, if y G [0, \) 

1, if y € [|,lj. 


9{y) = 


2.1. Numerical scheme. For the numerical scheme of equation (2.0.11 we use the 
following grid 

(xj,Vk) = (O' - 1 )A x,(k- 1 )A y),j = 1 , 2 ,--- ,M,k = 1,2,- •• , N, 

with Ax = £/(M — 1) and Ay = 1/(N — 1) the grid sizes in x— and y— direction, 
respectively. We denote the numerical approximation of u(xj, yk) by Uj^k- Similar 
to one-dimensional case, we use standard central differences for all derivatives in 


equation (|2.0.1|). The resulting nonlinear system of algebraic equations can be 

N{U) = AU + bb + qe u = 0, 


written in the form 
( 2 . 1 . 1 ) 
where 


( B 


A 


NxN = 


Ay 2 1 

o 


V o 


Ay 2 J 

B 


O 


O 

Ay 2 ^ 


Ay 2 1 

o 


B 

_2_ / 
Ay 2 J 


O \ 


o 

Ay 2 ^ 

B 
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5(M-2)x(M—2) - diag(-2 (• • • 


—2( 


1 

Ax 2 


1 

Ay 2 


), 0 ), 



(h\ 


( Ul ) 


( 0 ^ 


/ U2i \ 

bb = 

&2 

,u = 

U 2 

•> hi — 

0 

= 

U 3 i 


K b N) 


\UN ) 


\a ^givi)} 




Now we would 
suitable initial 


use Newton iteration scheme to solve equation (2.1.1). First give a 
guess. Solve 


J AU + bb + q = 0, 
l e u = l. 


we obtain 


U° = -A- 1 (bb + q-!t). 


The Newton iteration scheme to solve equation (2.1.1) is the following 


given, U°, 

solve J(U l )s l = —N(U l ), 
U l+1 = U l + s l , 


with 


U° 


—A 1 (bb + q ■ it), J(U) = A + q ■ diag{e u ), s l = 


N{U l ) 

J(U l ) ' 


2.2. Numerical results. Similar to numerical scheme for one-dimensional case, 
we use central difference method and Newton iteration method to solve equation 
(2.0.1). Let parameter q be 0.8, 0.5, 0.3 respectively, given 1 = 1, grid sizes 
Ax = Ay = 0.1, we could obtain the number of iterations, solutions and the 
residues as the following table and figures: 


q 

l 

Ax 

Ay 

number of iterations 

residue 

0.8 

1 

0.1 

0.1 

3 

1.37145172e-013 

0.5 

1 

0.1 

0.1 

2 

1.25775777e-009 

0.3 

1 

0.1 

0.1 

2 

2.03563927e-011 
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Figure 
1. q = 0.8 


Figure 
2. q = 0.5 



y 

Figure 


3. q = 0.3 


3. Conclusion 

We have considered a nonlinear boundary value problem describing the onset 
of a thermal explosion in a vessel in one-dimensional and two-dimensional forms. 
Numerical experiments show that Central different method combining Newton it¬ 
eration method is suitable to solve the thermal explosion model (0.0.21 if we find 
a suitable initial guess. In one-dimensional case (1.0.3), Richardson extrapolation 
shows that the convergence order of our algorithm is 2; in two-dimensional case 
(2.0.1), numerical experiments show that the number of iterations are less or equal 
to 3 and the residues are small. 
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